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APPLICATION OF BROUWER'S ARTIFICIAL -SATELLITE 

THEORY TO COMPUTATION OF THE 

STATE TRANSITION MATRIX 

By George H. Born and James C. Kirkpatrick 
Manned Spacecraft Center 

SUMMARY 


Brouwer's solution of the artificial- satellite problem without drag is used to 
obtain analytical state-transition-matrix expressions that include secular as well as 
long- and short-period effects of planetary oblateness. A comparison of the accuracy 
of several different models of the state transition matrix is made. These models 
include (1) a Keplerian model, (2) a model based on first- and second-order secular 
terms in Brouwer's theory, and (3) a model including first- and second-order secular 
terms in addition to first-order long- and short-period perturbations. It is demon- 
strated that the accuracy of the Keplerian model degenerates rapidly in comparison 
to either of the models based on Brouwer’s theory. In addition, it is shown with nu- 
merical results that including the effects of long- and short-period perturbations im- 
proves the accuracy of the transition matrix by one to two orders of magnitude over 
the secular model and three to four orders of magnitude over the Keplerian model. 


INTRODUCTION 


Determination of the state transition matrix is a necessity in orbit determina- 
tion and navigation theory. Numerous methods that can be used to compute the state 
transition matrix exist. These methods include direct integration of the variational 
equations, numerical differentiation, and various analytical formulations. For ex- 
ample, a Keplerian model for any type of conic motion has been developed by Good- 
year (ref. 1). Recently, Ditto (ref. 2) applied the Peano-Baker method to the 
integration of the variational equations. In this method, numerical integration is 
applied only to the departure from a simplified analytical model. 

The objective of this study was to obtain an analytical formulation of the state 
transition matrix that accounted for the effects of planetary oblateness for use in long- 
term satellite navigation studies. The method used was to differentiate the solutions' 
in Brouwer’s artificial-satellite theory (ref. 3) to obtain the transition matrix in 
terms of Keplerian elements. Although only the oblateness was considered in this 
study, the inclusion of long-period and secular effects for the remaining harmonics 
in Brouwer’s theory would be a simple extension. 



A numerical comparison of the accuracy of several different models of the 
transition matrix was made. These models included (1) a Keplerian model; (2) a 
model based on first- and second-order secular perturbations in £2, oo, and M and 
mean values of a, e, and I as given by Brouwer’s theory; and (3) a model includ- 
ing first- and second-order secular perturbations in £2, co, and M in addition to 
long- and short-period perturbations in all the orbital elements. The terminology 
first- and second-order perturbations means that the perturbations are proportional 
2 

to JgQ and JgQ , respectively. 

In this study, model 2 will be referred to as the secular model, and model 3 
will be called the complete model. Both models require the use of mean elements to 
be accurate through the first order. The reason epochal elements are not used in the 
secular model, which frequently is done erroneously, is discussed in the section on 
analysis. 

The authors gratefully acknowledge the assistance of Mr. Claude E. Hildebrand, 
Jr. , of The University of Texas at Austin, who helped with the analysis. 


SYMBOLS 


A matrix of time-dependent coefficients 

a semimajor axis 

a semimajor axis based on total energy 

B generic term for transition matrix 

e eccentricity 

second-order portion of doubly averaged Hamiltonian 
f true anomaly 

G (,a) 1/2 (! - e 2 ) V2 

H G cos I 

I • inclination 

JgQ oblateness coefficient 


2 



L (fia) 1 ^ 

M mean anomaly 

N dimension of the state vector 

n mean motion 

n mean of mean motion 

n mean motion based on total energy 

R relative error 

r radius 

r fi mean radius of planet 

S column vector of Cartesian coordinates 

t time 

V velocity 

X perturbation vector 

x, y, z Cartesian coordinates 
e column vector of orbital elements 

77 Vl - e 2 

e cos I 

ji gravitational constant 

<£>(t, t 0 ) state transition matrix 

£2 longitude of ascending node 

co argument of pericenter 

Subscripts: 
o initial value 

p evaluated at pericenter 

s secular rate 
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Operators: 


(• ) derivative with respect to time 

(' ) first-order averaged variables 

(”•) mean elements 

(*) evaluated on the reference orbit 

ANALYSIS 


Theoretical Development of State Transition Matrix 

The matrix differential equation for the state transition matrix is given by 

$(t,to) = A *(t,to) (1) 


where 



An element of $(t, t ) is given by 





( 2 ) 


(3) 


where (i, j = 1, 2, . . . , N), e i is a generic symbol for any of the orbital elements, 
and N is the dimension of the state vector. Also 


A = 


3e(t) 


(4) 


Equation (1) represents a system of N by N linear differential equations with known 
time-dependent coefficients. This so-called variational equations system may be 
solved by numerical integration. 
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One alternative to integration of the variational equations is the differentiation 
of an analytical solution, such as Brouwer's, with respect to the initial conditions to 
obtain the state transition matrix directly. Brouwer's theory separates the secular 
and periodic terms; therefore, contributions to the state transition matrix of the sec- 
ular and periodic portions also are computed separately. This separation permits the 
computation of a transition matrix based on any combination of secular and periodic 
terms desired. 

Because Brouwer' s theory is written in terms of mean elements and because 
differentiation with respect to epochal elements is desired, the chain rule must be 
used; that is, 


where 



0e(t) 

~0€" 

Ml 

r’tyj 

-M 


r-co) 


N‘o)i 

No). 




(5) 


( 6 ) 


and 


0a 8a 0a 

05^ 0e” ' • ‘ 0M" 



0M 0M 0M 

0a" 0e” 0M" 


The matrix 


0e 


"Co) 


3£ Co) 


is a constant matrix that needs only to be calculated once. E. 


the transition matrix in terms of Cartesian coordinates is desired, equation (5) must 
be multiplied by the appropriate transformation matrices as follows 
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where S represents the vector of Cartesian coordinates. The matrix 3S/3e is 
obtained by differentiation of the expressions that relate the Cartesian coordinates and 
the orbital elements. These expressions may be found in most fundamental celestial 
mechanics texts (e. g. , ref. 4). Because the transformation that relates the Cartesian 
coordinates to the orbital elements is unique, 


TO) 


TO 

TO 


TO 


(9) 


and this matrix is formed only once. 

The differentiation of Brouwer' s solutions with respect to the mean elements so 
that the results are accurate through the first order in the periodic terms and through 
the second order in the secular terms is straightforward but lengthy. Consequently, 
the results are not given in this study. However, a copy of the computer program 
containing the results will be supplied by the authors on request. 

The analytical models discussed are based on Keplerian elements and, there- 
fore, are unacceptable for circular orbits or for orbits with zero inclination. In addi- 
tion, Brouwer's theory is invalid within approximately ±2° of the critical inclination. 
The problems associated with zero eccentricity and inclination may be eliminated 
with a redefinition of the orbital elements as suggested by Lyddane (ref. 5). 


Computation of Mean Elements 

The mean elements for Brouwer' s theory are obtained accurate to the first 
order in JgQ by evaluation of Brouwer’s solutions at epoch and by replacement of 

mean elements with epochal elements in the first-order terms. For example, the 
mean value of a is given by 



Computing mean elements with epochal elements in the first-order terms results in 
second-order errors in the mean elements. As noted in Breakwell and Vagners 
(ref. 6), iteration of Brouwer's solution for the mean elements cannot reduce this 
error because the solutions do not contain second-order periodic terms. However, 
the secular portion of Brouwer’s theory is accurate to the second order; therefore, 
the mean value of a must be accurate to the second order to avoid second-order sec- 
ular errors being introduced in the mean value of the mean motion through the term 

M 1/2 ■ 

>2 
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The mean value of the mean motion correct to the second order may be computed 
according to Breakwell and Vagners (ref. 6) or Hildebrand (ref. 7). In terms of the 
total energy, ft is defined as 


1/2 

8 = ^75 HD 

<x 


where 


M _ M , 2 

2l" 2a~ “ 3 

o r 

o 


(- 








(12) 


Replacing equation (12) with its average through second order and using equa- 
tion (11) results in equation (13). 


n = n 


1 + 3 


m 2 k 2 


L'G* 




4 2 

,2 o M K„ 


1 3 H 


3L 



(13) 


The quantity Fg** is the second-order portion of the Hamiltonian averaged with 

respect to mean anomaly and argument of pericenter. It is given by equation (29) of 
Brouwer (ref. 3) as 


6 „ 2 
M K 2 

F ** = =— 

2 T ,10 


15 L* 



9H 


G' 


15 L’ 7 L 2H’ 2 
32 G’ 7 \ " G’ 2 



(14) 


Equation (1 3) now may be solved for n, which is the desired result. 

The use of mean elements is necessary to maintain the state transition matrix 
accurate to the first order. If epochal elements are used in a first-order secular 
model for the transition matrix, as commonly is done, a first-order secular error in 
the mean anomaly will exist. In fact, inclusion of the first-order secular term in the 
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computation of the mean motion may result in a transition matrix that is less accurate 
than a Keplerian model. To demonstrate this situation, the solution for the mean 
value of the mean motion through the first order is considered. 


1/2 

n- M 


1/2 2 t 

3 ^ r e J 20 


+ 4 


,,7/2 3 
a" r] 


(-1 + 30 2 ) 


To the first order 


a” 



(15) 


(16) 


Denoting the term of order J 9n by J ’ and substituting equation (16) into equa- 
tion (15) yields 


n = 


1/2 1/2 
( a o - J 2of /2 Vo - J 20’) 7/2 


V 


(17) 


Expansion of equation (17) and retention of first-order terms result in 


- 1/2 

n = ji 


-3/2 3 -5/2 T , -7/2 

a o + 2 a o J 20 +a o 


|^(. 1 + 3 „ 2 ) 

V 


°( J 20 2 ) 


(18) 


The first-order error introduced by the use of n = ^y^ is 


1/2 




1/2 j 3 -5/2 


2 % 


J 20’ + a o 


■7/2 


v 


(19) 


The first-order error introduced by use of equation (15) with the epochal value for a is 


^u 1 / 2 a' 5/2 J ’ 

2 M a o d 20 


( 20 ) 
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Therefore, if the initial conditions are such that the two error terms of equation (19) 
are of opposite sign, the Keplerian model will be more accurate than a secular model 
using epochal elements. 


RESULTS 


Discussion of Comparison Criteria 

The three analytical models considered, the Keplerian, the secular, and the 
complete model, were evaluated by comparison with numerical integration of the 
variational equations. The numerical integration was performed with a double- 
precision fourth-order Runge Kutta algorithm. The test cases used were a highly 
elliptical orbit around Mars and a near-circular orbit around the Earth. The Mars- 
orbital data were run because the results could be applied to long-term navigation 
studies for a Mars orbiter. All results presented are for transition matrices trans- 
formed to Cartesian coordinates. The disturbing force in all cases was planetary 
oblateness. 

It was assumed that Brouwer's theory for an Earth satellite is valid for a Mars 
orbiter. The fundamental assumption for Brouwer's theory was that all perturbing 
effects are second order compared to oblateness. 

To compare the accuracy of the various models, the relative errors as defined 
by Ditto (ref. 2) were used. The relative error is defined as follows. If 1IA|| is the 
norm of the transition matrix obtained by numerical integration of the variational 
equations; that is, 


HAII = £ £ IA I (21) 

i=l j=l 1 131 


and ||B|| is the norm of the transition matrix computed from any of the analytical 
solutions, then the relative error R is defined as 


R = 


II A - Bll 
"IfAll 


( 22 ) 



Discussion of Figures 

The relative error was computed at 0. 02-day increments for 25 days for the Mars 
satellite and at 10-minute increments for 14 400 minutes for the Earth satellite. The 
results for the Mars orbiter are shown in figure 1. The initial conditions for the Mars 
or biter are as follows. 


a = 0. 414245 x 10 8 feet 
o 

= 28.66° 

0 

e = 0. 7 

w =318.116 

0 

0 

I = 19.583° 
0 

o 

O 

II 

o 

2 


J 2Q = 0. 002011 


(23) 


The exponent of R is approximately equal to the number of digits of agreement 
between the analytical model and the integration of the variational equations. As seen 
from figure 1, the complete model agrees to four places; the secular model agrees to 
two places; and, after the first 2 days, the Keplerian model fails to agree even to one 
digit with integration of the variational equations. 

The relative errors for an Earth orbiter are presented in figure 2. The initial 
conditions for this orbit are as follows. 


a = 0.2378932 x 10 8 feet 
0 

II 

o 

a 

e = 0. 1 
0 


o 

O 

II 

O 

3 

I = 75° 
0 


M = 0° 
0 


J 2Q = 0.001083 

) 


(24) 


The agreement of the various models was approximately the same as for the Mars or- 
biter except that agreement to one more place was obtained with the secular model. 

The nominal trajectory for each model was generated with the state-propagation 
equations corresponding to the particular model. For example, the Keplerian results 
were based on a Keplerian nominal, while the secular-model nominal was generated 
from the following relationships. 


10 



a(t) 

= a” 


e(t) 

= e" 


I(t) 

= i" 


O(t) 

= ft" 

+ ft g At 

co(t) 

= co" 

+• co At 
s 

M(t) 

= M" 

+ M At 
s 


( 25 ) 


To demonstrate the magnitude of errors encountered by propagating a perturba- 
tion vector with the various models, a small perturbation vector was mapped with the 
analytical transition matrices and with the transition matrix obtained by numerical 
integration of the variational equations. These results were compared to actual inte- 
gration of the equations of motion with the perturbed initial conditions. 

For all results shown, the initial conditions were perturbed by the addition of the 
following perturbation vector to the initial state vector. 



500 

ft 

500 

ft 

500 

ft 

1.0 

ft/ sec 

1.0 

ft/sec 

1.0 

ft/ sec 


( 26 ) 


This perturbation vector was propagated for both the Mars and Earth satellite 
orbits. The results for the Mars orbit are shown in figures 3 to 16. The magnitudes 
of the position- and velocity -perturbation vectors for the Mars orbiter computed from 


AX(t) = *(t,t 0 )AX(t o ) 


(27) 


are shown in figures 3 and 10, respectively. The results shown in these two figures 
are computed from the Keplerian transition matrix; however, differences between all 
models are indiscernible when plotted on this scale. It should be noted that the magni- 
tudes of these perturbations grow very rapidly. Hence, the linear state transition ma- 
trix will propagate this perturbation vector accurately for only a short period of time. 
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The result shown in figure 4 is the difference in the magnitude of the perturbation 
vector as computed by the Keplerian transition matrices and the numerical integration 
of the equations of motion with perturbed initial conditions. The same result is shown 
in figures 5 to 7 for the secular model, the complete model, and integration of the var- 
iational equations, respectively. 

A comparison of figure 4 with figure 5 shows that the maximum error in 
the magnitude of the perturbation vector after one revolution (0. 5 day) is 10 nautical 
miles for the Keplerian model and 3 nautical miles for the secular model. The com- 
plete model and integration of the variational equations both have errors of about 
1. 2 nautical miles as shown by figures 6 and 7. The primary source of error for the 
Keplerian model is the error in mean motion, which results in a secular- error growth 
in the mean anomaly. Consequently, at a given time, the Keplerian transition matrix 
is computed at an incorrect orbital position. This problem is virtually eliminated by 
using a mean or average value for the mean motion as evidenced by the secular-model 
results in figure 5. 

The errors resulting from numerical integration of the variational equations are 
presented in figure 7. Neglecting the effects of numerical integration errors (round off 
and truncation), the errors shown in figure 7 are caused solely by nonlinearities. It 
should be noted that the error histories for the complete model and for the variational 
equations are identical because the transition matrices agree to four places as shown 
by figure 1. 

A plot of the magnitude of the angle between the perturbation vectors in position 
computed from the Keplerian model and from integration of the equations of motion is 
presented in figure 8. The corresponding results for the variational equations are pre- 
sented in figure 9. The results for the complete model were identical to those of the 
variational equations. Because the angular deviations for the secular model were only 
slightly larger than those for the complete model, the deviations are not shown. Once 
again, the primary source of error for the Keplerian model is the use of an improper 
mean motion. 

Figures 10 to 16 correspond to figures 3 to 9, respectively, except that these 
results are presented for the velocity -perturbation vector. The same comments made 
for the position-perturbation-vector errors apply to these figures. 

The magnitude of the position-perturbation vector as propagated by the Keplerian 
transition matrix for the Earth orbiter is shown in figure 17. The errors in the mag- 
nitude of this vector for the Keplerian and complete models are given in figures 18 and 
19, respectively. The results from the complete model were virtually identical to the 
results from integration of the variational equations. Errors in the secular model were 
larger than for the complete model but much smaller than those of the Keplerian model. 


Theoretical Analysis of Results 

It should be noted from figures 3 and 10 that the magnitude of the position- 
perturbation vector has grown to approximately 190 nautical miles in position and 
800 ft/ sec in velocity after 0. 5 day, which is the pericenter passage of the first revo- 
lution. Although seemingly rather large, these numbers are results of the large 
eccentricity and can be verified easily analytically. Assume that it is desired to 
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compute the perturbations in position and velocity at successive pericenter points 
because of a small perturbation in the state vector at the initial time. The pertur- 
bation in the radius vector Ar p will have components Ar and r Af along and 

normal to the radius vector, respectively. At the pericenter, r = a(l - e); therefore, 

0f 0f 

Ar = Aa (1 - e) - a Ae. The component r Af is now computed from df = -^ de + dM 
where 

H=(777 + ?) slnf < 2 «> 

and 


af 

3M 



(29) 


Therefore, for small perturbations in the orbital elements, 

Af = ( — I—* + sin f Ae + ^ \fl - e 2 AM 
\1 - e / r 

2 + e cos f . f . 

1 - e 

+ (l - e 2 ) ^ (1 + e cos f) 2 (aM q - Aa t) (30) 

When r Af is evaluated at pericenter, 

/i pN 1 / 2 

(r Af) p = (a AM q - 377m Aa) (m = 0, 1, 2, . . . ) (31) 


where m denotes the number of pericenter passages. Because of the dominance of the 

secular term in (r Af) , Ar may be written as follows. 

P P 
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3irm Aa 


(32) 


Ar p 



1/2 


An expression for Aa in terms of the perturbations in r and V can be obtained 

V 2 

from the energy expression " 2 a = “ 2 “ " Equation (33) results when the first varia- 
tion is determined. 


2a 2 V AV 2a 2 

+ ~T 
r 


Ar 


(33) 


where 


V AV = xAx+yAy + zAz 


(34) 


and 


Ar = p (x Ax + y Ay + z Az) 


(35) 


An expression for 

ating for Ar . 

P 



is obtained by solving equation (33) for 


AV and evalu- 


AV = 
av p 2V 


2 Ar 


(36) 


When a substitution is made for 


Ar^, the following equation results. 



JL- 

2V 


Aa . 

f 1+6 l 

l^ 2 6vm Aa 

La 2 

U - e) 

1 2 
r 


( 37 ) 
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When the term 



is neglected, the following equation results. 


4V p 


fl + e\^ 3jrai Aa 
U - e) r 2 v 


(38) 


The magnitudes of the position- and velocity-perturbation vectors at pericenter com- 
puted with equations (32) and (38) are 185 nautical miles and 775 ft/sec, respectively, 
which are in very good agreement with the results shown in figures 3 and 10. 

CONCLUDING REMARKS 


It has been shown by a numerical example that Brouwer's artificial- satellite 
theory can be used to generate transition matrices that agree closely with numerical 
integration of the variational equations for long periods of time. It has been demon- 
strated that inaccuracies in the mean motion cause the accuracy of a Kepler ian transi- 
tion matrix to degenerate rapidly. However, the accuracy of a secular model based on 
mean elements or of a complete model containing secular and periodic terms did not 
decrease for the time period considered. Although the complete model, of course, is 
more accurate than the secular model, using .the complete model involves many more 
terms and requires more computer time. Hence, the secular model should be used 
whenever possible. 


Manned Spacecraft Center 

National Aeronautics and Space Administration 
Houston, Texas, May 27, 1970 
914-50-17-08-72 
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Figure 3. - Position-perturbation magnitude for Keplerian model, 

Mars orbiter. 
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Figure 8. - Position-perturbation angular error for Keplerian model, 

Mars orbiter. 
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Figure 15. - Velocity-perturbation angular error for Keplerian model, 

Mars orbiter. 
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Figure 19. - Position-perturbation error for complete model, 

Earth orbiter. 
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